Thermodynamics of ultra-small metallic grains in the auxiliary-field Monte Carlo 
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We use an auxiliary-field Monte Carlo (AFMC) method to calculate thermodynamic properties 
(spin susceptibility and heat capacity) of ultra-small metallic grains in the presence of pairing cor- 
relations. This method allows us to study the crossover from bulk systems, where mean-field BCS 
theory is valid, to the fluctuation-dominated regime of ultra-small particles at flnite temperature. 
The computational effort at low temperatures is significantly reduced by exploiting a simple renor- 
malization method. 
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I. INTRODUCTION 

The properties of ultra-small metaUic grains (nanopar- 
ticles) differ from those of bulk superconductors. Larger 
grains are well described by the Bardeen-Cooper-Scrieffer 
(BCS) mean- field theory^. However, as the size of the 
grain becomes smaller, fluctuations around the BCS 
mean- field solution become important. This crossover 
from the BCS limit to the fluctuation-dominated regime 
is controlled by the ratio A/(5 between the bulk BCS gap 
A and the mean level spacing i5. BCS theory is valid for 
A/5 ^ 1, while fluctuations become important as this 
ratio decreases and becomes similar to or smaller than 

Most of the recent interest in such ultra-small metal- 
lic grains has been motivated by experiments in the 
mid 1990s^ in which the spectra of individual aluminum 
grains were determined by measuring the non-linear con- 
ductance. These spectra were found to depend on the 
magnetic field strength as well as the number parity 
of electrons in the nanoparticle. Pairing correlations 
were clearly observed in the larger metallic particles, and 
the results were explained qualitatively using a number- 
parity projected BCS model.— However, for smaller par- 
ticles with nominal radii in the range ~ 2 — 5 nm it was 
not possible to distinguish between the spacing due to the 
discreteness of the single-particle spectrum and the pair- 
ing gap. The authors of Ref. |3i emphasized that this does 
not mean that pairing correlations cannot be observed 
in these samples but rather that conductance measure- 
ments might not be the appropriate tool to study the 
fluctuation-dominated regime. 

It was suggested that number-parity effects of pairing 
correlations in the fluctuation-dominated regime can be 
observed in thermodynamic quantities such as the spin 
susceptibility;^, Using a combination of different meth- 
ods, the authors of Ref. 5 found a reentrant behavior 
of the spin susceptibility of grains with an odd num- 
ber of electrons. For low temperatures they exploited 
an exact solution of the BCS pairing Hamiltonian known 
as the Richardson solutioni^ Results were presented for 
A/5 < 1 and T < 1.6 5 because the Richardson solution 



becomes less tractable at larger values of A/5 and/or 
higher temperatures. Their results at high temperatures 
were based on the static-path approximation^ which ig- 
nores quantum fluctuations. Signatures of pairing cor- 
relations in the spin susceptibility of ultra-small grains 
were also studied in Ref. H. 

Here we present a method to study thermal quantities 
such as spin susceptibility and heat capacity in the full 
temperature range and for arbitrary values of the ratio 
A/5. We use the auxihary-field Monte Carlo (AFMC) 
method developed in the framework of the interacting 
nuclear shell model, which is known (in nuclear physics) 
as the shell model Monte Carlo (SMMC) method.^-^o The 
SMMC approach was used to study thermal signatures of 
pairing correlations in the heat capacityiiii^ and the mo- 
ment of inertia^'^ of nuclei. Similar methods were used in 
the study of strongly correlated electron systemsJ^ Here 
we adapt the AFMC method to study pairing correlations 
in metallic nanoparticles. The dimension of the many- 
particle space depends combinatorially on the number of 
single-particle levels and/or number of electrons, and a 
direct diagonalization of the many-body Hamiltonian be- 
comes impractical. AFMC scales as a power law in the 
number of single-particle orbitals and can be carried out 
in very large model spaces. Recently, a different Monte 
Carlo approach that is based on canonical non-local loop 
update algorithm and is specific for a pairing Hamilto- 
nian was used to study metallic nanoparticles 

The thermodynamic properties of grains with chaotic 
single-particle dynamics are expected to display meso- 
scopic fluctuations because of the randomicity associ- 
ated with their single-particle Hamiltonian. For example, 
mesoscopic fluctuations in the orbital magnetization of a 
chaotic grain in the presence of an orbital magnetic field 
were studied in Ref. [H. Our studies here are focused 
on grains with equally-spaced spectra and do not include 
the effects of mesoscopic fluctuations. 

The outline of this paper is as follows. In Sec. [H] we dis- 
cuss the AFMC method and its particular realization for 
the BCS pairing Hamiltonian. The method is based on 
a Hubbard-Stratonovich transformation,^^ and the num- 
ber of particles is fixed using a particle- number projec- 
tion. Section IIIII discusses various ways to renormalize 
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the pairing coupling constant when the band width is re- 
duced. By reducing the band width, the low-temperature 
AFMC calculations can be carried out more efficiently. In 
Sec. IIVI we use AFMC to calculate the spin susceptibility 
and heat capacity of grains with different values ol A/8. 

II. AUXILIARY-FIELD MONTE CARLO 
METHOD 

We briefly discuss the model and then describe the ap- 
plication of the AFMC method to an ultra-small metallic 
grain. 



electrons moving in external time-dependent auxiliary 
fields. This is accomplished by dividing the time in- 
terval (0, (3) into Nt time slices of length A/3 each 

such that e^^^ — {e~^^^^ *. For each time slice 

-A/3H ^ -A|3Y.i<^^fHX[..p-^Pi[{P^^+P^j? -{P^j-Pij?] tO 

order (A/3)^. The interaction terms in the exponent ap- 
pear as sum of squares and the propagator for one time 
slice can be written as Gaussian integrals over real aux- 
iliary fields and crf^ 

g-A/5f [(p.,+p,,)^ -(P..~P,,r] f daedal (3) 

A-K J J J 



A. Model 

We describe the metallic grain by the BCS pairing 
Hamiltonian^. 



i = -Na 



No 



-g 4+4-Cj-cj+ A'^) 

l,j = -Na 



where cl^ is the creation operator for an electron in the 
single-particle level i with either spin up (tr = or spin 
down (cr = — ). We assume 2No + 1 levels with a cutoff 
No determined by the Debye frequency uju via No ~ 
ujd/8 where 6 the mean single-particle level spacing. The 
second term in ([T]) is a pairing interaction, which scatters 
only time-reversed pairs of electrons. In the following 
we will assume an equidistant spectrum e, = iS {i — 
—No, . . . , No) and half filling, i.e., the number of electrons 
is either N = 2No (even) or TV = 2No + 1 (odd). The 
coupling constant g is related to the BCS pairing gap 
(at zero temperature) through A = 2lli-o^^^^^ ■ Both the 
single-particle mean level spacing and the pairing gap are 
assumed to be small compared with the Debye frequency, 
i5, A <tJD- 



B. Hubbard-Stratonovich transformation 

The BCS Hamiltonian ([1]) can be written in a density 
decomposition as 
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X^KPy + P^3f - (P^j - Pu)^] ' (2) 



where fii = cl^Ci^ + c|_Ci_ is the occupation number 
operator of level i, pij = c\j^Cjj^ are spin-up density op- 
erators and pij = c\_Cj- is the time reverse operator of 
Pij describing a spin-down density operator. 

The operator e^^^ at inverse temperature (3 — 1/T 
can be viewed as a many-body propagator in imag- 
inary time [3. In the Hubbard-Stratonovich (HS) 
transformation it can be written as a path integral 
over one-body propagators that describe non-interacting 



Using different complex fields at each time slice 
<^i]{Tn) = o-^jiTn) + jcr,fj(T„) (r„ = nA/3 with n = 
I, . . . , Nt), the propagator e~^^ can be written as a path 
integral 



e-^"^ I V[a]G{a)U,{(3,0) 
with the measure 



v[a] = n 

and Gaussian factor 
G((t) = exp 



daij{Tn)da*j{Tn) Af3g 



2i 



47r 



(4) 



(5) 



-f A;3^|a.,(r„) 



(6) 



Ua in ll]) describes the propagator of non-interacting 
electrons in external time-dependent auxiliary fields 
<T^J (r) 

{f3^ 0) = e-'^'^'*"(^"*) . . . e-^'^'^'f^i) (7) 
with a time-dependent one-body Hamiltonian 

'^'^(t) = ~ 2 XI '^^3 + ('^)] • 

i ij 

Eq. ([U is known as a Hubbard-Stratonovich (HS) 
transformation. 

Another HS transformation can be obtained using 
good-spin density operators 



Pkm 



(U) = E ( V2, a./2; 1/2, a,/2\KM) c^5,.^. (9) 



with c,a, = (-I)(^+'"')/2ci_^^ and K = {),!; -K < M < 
K . We rewrite the pairing Hamiltonian ([T]) in terms of 
these good-spin density operators as 
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where pkm = i^'^) Pk-m is the time-reversed oper- 
ator of Pkm- Eq. pop can be rewritten in a form similar 
to (121) but now with pKM{ii) replacing pij. An HS de- 
composition can be carried out, introducing the complex 
fields <7KM{ij)- This HS decomposition is similar to the 
one used in SMMC. 

To account for finite-size effects in the grain, it is neces- 
sary to evaluate thermal expectation values in the canon- 
ical ensemble with fixed electron number N . For an ob- 
servable O 
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(11) 



where Tr^v denotes a trace at a fixed number 
of particles N. Here = G(cr)|TrArC/^(/3, 0)| 

is a positive definite weight function, — 
Trjv/7<,(/3,0)/|TrArt7<,(/3,0)| is the "sign" function, 
and {0)a = TiN[OUa{f3,0)]/TrNUa{P,0). Both Tr^C/, 
and (0)cr can be evaluated using matrix algebra in the 
single-particle space. For example, the grand-canonical 
trace of the one-body propagator C/(j(/3, 0) in Fock space 
is given by 



Tr{/,(/3,0) =det[l-|-U,(/3,0)] , 



(12) 



where U<^(/3, 0) is the (4iVo + 2)x (ANo + 2) matrix rep- 
resenting Ua{P,0) in the single-particle space (each level 
is spin degenerate and thus the dimension of the single- 
particle space is twice the number of levels) . 



C. Particle- number projection 



D. Monte Carlo sign 

The functional integral in (jlip is a multi-dimensional 
integral over a large number of auxiliary fields. In the 
AFMC method, this integration is carried out by Monte 
Carlo methods. The auxiliary fields are sampled accord- 
ing to the distribution Wa, and for AI samples {fXi}, the 
thermal expectation value (jlip is estimated by 



(16) 



In general, the sign function (f>cr is not necessarily pos- 
itive and fluctuates from sample to sample. When the 
statistical error of Re$o- is larger than its average value, 
the observables cannot be reliably estimated, a problem 
known as the Monte Carlo sign problem. In the case 
of the attractive pairing interaction, the single-particle 
Hamiltonian ((8]) is time-reversal invariant, ha- = ha-, and 
so is Ua- Therefore the eigenvalues of Ucr(/3,0) come in 
complex conjugate pairs (corresponding to an eigenstate 
and its time- reversed state). Thus the grand-canonical 
many-particle partition TrJ7cr(/3, 0) in is positive def- 
inite. When projected on an even number of electrons 
TrjvC/cr(/3, 0) remains almost always positive and there is 
no sign problem. When projected on an odd number of 
electrons TTNUa{P,0) can be negative and a sign prob- 
lem occurs at large values of (3. In practice, we use the 
reprojection methodi^ We carry out the Monte Carlo 
sampling with an even number of electrons and then re- 
project on an odd number of electrons. 



The canonical trace can be evaluated from TinX ~ 
Ti{PnX) where Pjy — d{N — N) is the particle-number 
projector. Since the number operator N takes a finite 
number of discrete values 0, 1, • • • , iNo 2, Pjv can be 
written as a discrete Fourier transform 
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4No + 2 
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(13) 



where = 2TTm/ {ANo -I- 2) are quadrature points. The 
canonical trace of Ua is then given 



(14) 



xdet[l + e''^™U^(/3,0)] . 
Similarly for a one-body observable 

4Af„+2 



Tr 



N 



(15) 



det[l-he"^'"U<,(/3,0)]tr 



l + e-»*'"U^i(/3,0) 



O 



where O is the matrix representing the operator O in the 
single-particle space. 



E. Discretization effects 

The discretization of imaginary time introduces a sys- 
tematic error in the evaluation of the HS integral 
We corrected for these discretization effects by carrying 
out AFMC calculations for several values of A/? (usually 
1/32 and 1/64 5^^) and then performing a linear extrap- 
olation to A/? = 0. 



III. RENORMALIZATION 

The AFMC method has the advantage that it only re- 
quires matrix algebra in single-particle space. The com- 
putational time scales as a power law ^ with the 
number of single-particle states and not exponentially as 
is the case with conventional diagonalization methods. 
We can further reduce the computational effort at low 
temperatures by using the Hamiltonian ([T]) in a trun- 
cated model space with a smaller cutoff Nr < No and a 
renormalized coupling constant g^- For a renormalization 
group (RG) approach that describes the BCS instability 
in the bulk, see Ref. 21. In this section we discuss various 
ways to determine gr for a finite-size system. 
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In the bulk, thermodynamic quantities are expected 
to be universal functions of T/A where A is the BCS 
gap. The T = BCS gap is determined from 1/g — 
l/{2Ek) where Ek = \J {et — A*)^ + A^ are the quasi- 
particle energies. Taking the bulk limit ^ — > of a 
picketfence spectrum in a band between — and lou 
with half filling (/x w 0), the solution for A is given by 
A = u;£)/sinh(l/A) where A = gjb, or A = 2ixioe~^l^ 
for A <C 1. An effective Hamiltonian in a truncated 
flat band of half width D (smaller than lod) can be 
obtained by rcnormalizing the coupling constant A to a 
value Ar such that the gap parameter remains fixed, i.e., 
A = L>/sinh(l/A^). 

In a finite-size system, the discreteness of the single- 
particle spectrum becomes important and there are two 
energy scales: the BCS gap A and the single-particle 
mean level spacing 5. When truncating to a band with 
2Nj. -I- 1 levels, we renormalize g such that (at fixed S) 
the BCS gap A, determined by the BCS equation for a 
discrete spectrum, remains fixed. Using the gap equation 
for both untruncated and truncated spectra, we find the 
following relation between the bare coupling constant g 
and the renormalized constant g^ 



f 
9 



1 



E 



(17) 



In the sum within the shell iV^ < < iVo (outside the 
truncated band), we have ignored the gap A assuming 
A <C Nr^. We have also assumed that the chemical po- 
tential is the same for No levels and for levels. For 
the picketfence spectrum at half filling this is exact by 
symmetry; /i = —(5/2 for an even number of particles TV 
and = for odd N . 

Equation P7)l can be rewritten as 



5 1-5 
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N^<\k\<N„ 
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2|e, 



Ml 



(18) 



This equation was also obtained using a perturbative 
treatment^ by taking the limit and wd — >• cxd 

at fixed A and b. A quite accurate estimate of the dis- 
crete sum in (|18|) for a picketfence spectrum is given 

by (for = 0) EAr,<fe<Ar„ Vl^l « InIHiI'^^I'' = 
ln[(2iVo + 1)/(2A^^ 1)]. Using this relation in (HH]) and 
A = 2ujD&^^l^ with UJD = {No + l/2)(5, we arrive at 
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2Af, 



A/5 



(19) 



For a picketfence spectrum, rather than using (jlSp . 
we can calculate quite accurately the discrete sum in 
the BCS gap equation using I/5 = ^k<\N \ l/(2^'fc) ~ 



(A/5)2. We find 
1 

arcsinh (^^f^ 



(20) 




FIG. 1: Renormalized coupling constant g,, versus A/ 5 for 
Nr = 5 (upper curves) and A^r = 10 (lower curves). Solid 
circles correspond to keeping the discrete BCS gap A fixed 
(at fixed S) and open circles describe the renormalization of 
Eq. p8|l . Solid lines correspond to ((20} and dashed lines are 
the bulk renormalization (|19|l . The inset shows the region 
A/5 < 0.5. 



Figure [T] compares the various renormalized values of 
the coupling constant gr versus A /S for truncated bands 
of Nr = 5 (upper curves) and Nr — 10 (lower curves). 
Open circles describe the renormalization of Eq. (|18p . 
while solid circles correspond to keeping the discrete BCS 
gap fixed. Dashed lines are the bulk renormalization (fTO)) 
and solid lines describe the normalization given by (j20p . 
We observe that for A/S > 0.5, the renormalization ([20]) 
essentially coincides with keeping the discrete BCS gap 
fixed, while the bulk renormalization (jl9p coincides with 
the renormalization (|18p of Ref. [13. The renormalization 
(fT8)) deviates from the discrete BCS renormalization at 
large values of A/S but this deviation becomes smaller 
as Nr increases. This is consistent with the condition 
A <^ NrS. A different behavior is observed for A/S < 0.5 
(see inset of Fig. [1]) where the renormalization ([20]) co- 
incides with ([TS]) . To test the different renormalization 
methods, we show in Fig. [5] the excitation gaps SEs for 
spins S = 1/2,3/2 and 5/2 (for an odd number of par- 
ticles) as a function of the width Nr of the truncated 
band for a fixed value of A/S = 3. These spin gaps were 
determined by solving Richardson's^uations. Open cir- 
cles correspond to Eq. (|18p of Ref. 22 while solid circles 
correspond to keeping the discrete BCS gap fixed. We 
observe faster convergence for the second method. How- 
ever, the discrete BCS equation is meaningful as long as 
A/S > 0.5. For A/S < 0.5, the method of Ref. (or al- 
ternatively Eq. ([^ for a picketfence spectrum) should be 
used and it gives the correct limit gr ^ when A/(5 — > 0. 

Figure [2] also demonstrates that the T = renormal- 
ization of the coupling constant approximately preserves 
the low-energy excitation spectrum of the grain. Thus, 
the same renormalization is expected to work well at fi- 
nite temperature T, as long as Nr is sufficiently large. 
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FIG. 2: The excitation gaps SEs (i.e., lowest excitation en- 
ergy for a given spin S) for spins S — 1/2, /3/2, 5/2 versus 
the half width A^r of the truncated band. The gaps are shown 
for A/ 5 = 3. Solid circles correspond to keeping the discrete 
BCS gap fixed while open circles correspond to the renormal- 
ization'^^ of Eq. (|18p . Dotted lines describe the asymptotic 
values of SEs for large Nr. 



At higher temperatures, higher-energy configurations be- 
come populated and it is necessary to increase the band 
width Nr- In our calculations we ensure that our model 
space is large enough by estimating the minimal required 
value of Nr at each temperature from the canonical Fermi 
gas results. We also note that the renormalization pro- 
cedure works well without introducing new coupling con- 
stants beyond g. For grains with equally-spaced spectra 
and given A/S, the various thermodynamic properties are 
expected to be universal functions of T/S, irrespective of 
the bandwidth. 

As an example we consider an ultra-small aluminum 
grain with A/S — 1, where A denotes the bulk pairing 
gap at zero temperature. The measured Debye frequency 
is « 34meV. Using the experimental value for the 
zero-temperature gap of thin films A « 0.38 meV, we can 
estimate the bare cutoff to be No « 89. The bare cou- 
pling constant is then given by g/6 ~ 1/ ln(2ti;£)/A) — 
0.193. In the actual calculations, we used truncated val- 
ues Nr — 10, 15, 25 (depending on temperature) and used 
the appropriate renormalized values gr of the coupling 
constant. 



IV. THERMODYNAMIC PROPERTIES 



We used AFMC to calculate thermodynamic proper- 
ties of a metallic grain in the crossover regime. In par- 
ticular, we calculated the spin susceptibility and heat 
capacity for grains with both even and odd number of 
electrons. 




FIG. 3; Spin susceptibility x (in units of the Pauli suscepti- 
bility Xp = 2/is/5) versus T/S for A/S = 1. Solid circles are 
the AFMC results for an even number of electrons and open 
circles are AFMC results for an odd number of electrons. The 
solid lines are obtained by using Richardson solution for all 
energy levels below an excitation energy of ~ 305. For com- 
parison we also show the spin susceptibility of the free Fermi 
gas in the canonical ensemble (dottes lines), and of the BCS 
mean-field theory (dashed line). 



A. Spin Susceptibility 

The spin susceptibility is defined as the magnetic re- 
sponse of the grain to an external Zeeman field H. Tak- 
ing the z axis along the field direction, the field couples 
directly to the spin projection Sz, and the spin suscepti- 
bility is given by 



= {{Si) - {Szf 



(21) 



where M is the magnetization of the grain, /is is the 
Bohr magneton and we have used a g-factor of 2 for 
free electrons. The spin projection is given by Sz — 

AFMC results for the spin susceptibility of a grain with 
A/5 = 1 are shown in Fig. [3l Solid circles (with statis- 
tical errors) are for even number of electrons and open 
circles are for odd number of electrons (at half filling). 
The spin susceptibility is measured in units of the Pauli 
susceptibility Xp — 2^1/^ (see below). We used Nr = 10 
and M = 16000 samples for T < 0.8 (5 (at T = 0.4 ^ we 
have used M = 64000 samples), Nr = I'o and M = 4000 
samples for 0.8^ < T < 1.5 5, Nr = 25 and M = 4000 
at higher temperatures. The larger number of samples 
at lower temperatures is necessary to overcome a mild 
sign problem for the reprojection on the odd number of 
electrons. We also solved Richardson's equations using 
the method of Ref. up to an excitation energy of 30 5, 
and the results are shown by the solid lines in Fig. O 
They agree with the AFMC results up to T « 1.25 5. At 
higher temperatures, it is necessary to increase the max- 




imal excitation energy. However, this becomes a difficult 
problem as the number of many-body levels proliferates 
combinatorially with excitation energy.^ 

For comparison, we have also calculated the spin sus- 
ceptibility of a free Fermi gas in the canonical ensemble. 
The results are shown by the dotted lines in Fig. [3] for 
even and odd number of electrons. 

At high temperatures, we can calculate x for & free 
Fermi gas in the grand-canonical ensemble. We have 

x(T) = 2/3mIE/^(1-/') (22) 

i 

where fi = 1/ [l -I- e'^''^'~'''] is the Fermi-Dirac occupa- 
tion of level i. For the half- filled picket-fence spectrum, 
/i = 0, and for temperatures T ^ 5 (but much smaller 
than the Fermi energy) 

This value of x is the Pauli susceptibility Xv We observe 
that the canonical Fermi gas susceptibility (dotted lines 
in Fig. [3]) essentially coincides with (l23l) (i.e., x ~ Xp) 
already for T/6 > 1. The AFMC results (that includes 
the pairing interaction) approaches Pauli's limit at high 
temperatures. 

In the BCS limit, we can calculate x from ([22]) but 
now fi are the quasi-particle occupation numbers fi = 
1/(1 -He''^') with E, = ^(e, - ^)2 -t- A2 being the 
quasi-particle energies. The BCS spin susceptibility is 
shown by the dashed line in Fig. [S] We observe that 
pairing correlations in x persist up to temperatures that 
are higher than the BCS critical temperature. 

For larger gap values, even-odd effects extend to higher 
temperatures and it is difficult to use Richardson method. 
On the other hand, AFMC calculations remain tractable. 
The AFMC spin susceptibility for A/S = 3 is shown 
in Fig. m Here we used Nr = 10 and M = 64000 



samples for T < 1.25 5 (at T = 0.67(5 we have used 
M = 384000 samples), Nr = 15 and M = 16000 sam- 
ples for 1.25(5 <T< 1.6 S, N,. = 25 and M = 4000 at 
higher temperatures. 

We observe the following: 

• The pairing interaction suppresses the spin sus- 
ceptibility when compared with the canonical free 
fermi gas susceptibility for both even and odd num- 
ber of electrons. At higher temperatures, the re- 
sults for the even and odd number of particles co- 
incide and approach asymptotically the Pauli spin 
susceptibility. 

• The spin susceptibility for an odd number of elec- 
trons shows a characteristic reentrant behavior, in 
agreement with findings in Refs. 5 and 23. The 
unpaired electron leads to a Curie-like divergence 
X ~ 1/T in the limit T ^ (even for the canonical 
Fermi gas). However, pairing correlations suppress 
the spin susceptibility, leading to a minimum in x 
at a certain temperature. This behavior of x is sim- 
ilar to the the behavior observed in the moment of 
inertia of an odd-even or odd-odd nucleus. ^'^ 

• In the presence of pairing correlations, number- 
parity (i.e., even-odd) effects already appear at 
temperatures T < 1.5(5 for A/S — 1, as compared 
to T < 1.0 (5 for the free Fermi gas. For larger values 
of A/S, the region of reentrant behavior is shifted 
to higher temperatures and is more pronounced. 
The spin susceptibility for an even number of elec- 
trons gets closer to its BCS approximation at larger 
values of A/S. 

B. Heat Capacity 

The heat capacity C = dE/dT is calculated in AFMC 
as a numerical derivative of the thermal energy E{[3) = 
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FIG. 6; Heat capacity for A/(5 = 3. See Fig. |3]for notation. 



Figs. [5] and [6] The BCS heat capacity is calculated from 
C ~ TdS/dT using the quasi-particle occupations for fi 
in ((25|l and is shown by the dashed lines. 
We observe the following: 

• For intermediate temperatures, the heat capacity is 
enhanced for an even number of particles as com- 
pared with the heat capacity for an odd number 
of particles. This effect is clearly seen for A/i5 = 3 
where the heat capacity for an even number of elec- 
trons displays a shoulder. A similar behavior was 
found in the heat capacity of even-even neutron- 
rich nuclei.— This even-odd effect is a signature 
of pairing correlations in the crossover regime. At 
very low temperatures this behavior is reversed and 
the heat capacity is more strongly suppressed in the 
even case (see the solid lines in Fig. [5|). 



Tr(iJe-'^^)/Tre-'3^. We use 

C^~P^^^^±^^t_EiP^^Om^. (24) 

The statistical errors are reduced by using the same set 
of auxiliary fields a for the calculation of both E{/3zt 5(3) 
and taking into account correlated errors.-^ Figures [5] 
and [6] show the AFMC heat capacity for A/5 = 1 and 
A/5 = 3, respectively, for both even (sohd circles) and 
odd (open circles) number of electrons. The sohd lines 
in Fig. [5] are Richardson results when energy levels up 
to '--^ 30 (5 are included. They agreed with AFMC for 
T < 1.3(5 but deviate at higher temperatures because of 
truncation effects. 

For a non-interacting Fermi gas in the grand-canonical 
ensemble, we calculate the heat capacity from C = 
TdS/dT, where the entropy S is given by 

5 = -2^[/,ln/,-H(l-/,)ln(l-/,)] . (25) 

i 

The factor of 2 in ([25|l accounts for spin degeneracy. In 
the high-temperature limit T ^ S and for a picketfence 
spectrum, we find 

The canonical Fermi gas heat capacity (for an even 
number of electrons) is shown by the dotted lines in 



• The AFMC heat capacity is suppressed with re- 
spect to the BCS heat capacity. For A/6 = 3, the 
shoulder in the heat capacity for the even number of 
electrons occurs around the BCS critical tempera- 
ture. However, for A/6 = 1, this shoulder structure 
occurs at higher temperatures than the BCS criti- 
cal temperature. This suggests that the even-odd 
effects occur on the scale that is the larger between 
A and 6. 



V. CONCLUSION 

We used an auxiliary-field Monte Carlo approach 
to calculate thermal observables of ultra-small metal- 
lic grains. The method is particularly useful for grains 
in the crossover regime between the BCS limit and the 
fluctuation-dominated limit, where the BCS gap is com- 
parable to the single-particle mean-level spacing. For an 
attractive pairing interaction, there is no Monte Carlo 
sign problem for an even number of electrons, and accu- 
rate calculations are possible. The computational effort 
can be further reduced by renormalizing the interaction 
in a truncated band. For a picketfence spectrum and for 
a given ratio A/ 6 of the pairing gap to the mean level 
spacing, the spin susceptibility and heat capacity are uni- 
versal functions of T/6. 
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